### R code to accompany Andean community model
### .csv files are available at Data Dryad Repository

rm(list=ls())
require(dclone)
require(R2jags)
library('rjags')
library('runjags')
library('doParallel')
library('random')
n.cores <- 3
setwd("C:/UserDirectory")

### covariate data
covar.data<-read.csv("Andes.covariate.csv") # covariate data for each site 
names(covar.data)

### observer covariates
OBSVAR<-as.factor(covar.data$observer) 
height<-covar.data$canopy.std
time<-covar.data$min.std
time2<-time*time

# site covariates
MF<-covar.data$MF
SG<-covar.data$SG
SC<-covar.data$SC
SUN<-covar.data$SUN
FFRS<-covar.data$FFRS
LFP<-covar.data$LFP
PAR<-covar.data$PAR
elev<-covar.data$elev.ST
elev2<-elev^2
lat<-covar.data$lat.std
tree<-covar.data$treecoverST
precip<-covar.data$Wprecip.meanST
seas<-covar.data$Wprecip.seasST
cloud<-covar.data$Wcloud.meanST

VAR<-as.matrix(cbind(SG,SC,SUN,FFRS,LFP,PAR,elev,elev2,tree,precip,seas,cloud,lat)) # linear relationships except for elevation
nVAR=dim(VAR)[2]  # number of continuous covariates on abundance

### species matrix
y<-read.csv("Andes.species.csv") # total detections by site and species
y<-as.matrix(y[,2:24])
nind=sum(y) # total number of observations
nsp=dim(y)[2] # number of species in data set
nSeg=dim(y)[1] # number of transects
Nin=y+1 # initial values for N

### dclass
spp.data<-read.csv("Andes.dclass.csv") # distance class file for each observation
dclass<-spp.data$dclass	# distance class for all observations
species<-spp.data$spec # species index for all observations
site<-spp.data$site # site index for all observations
db<-seq(0, 40, 5) # breaks for distance categories 
nG=length(db)-1 # number of distance categories
pix<-c(0.01563, 0.04688, 0.07813, 0.10938, 0.14063, 0.17188, 0.20313, 0.23438) # %area covered by each category
v=5 # width of distance categories

# new variable for elevation covariate predictions
m<-seq(-2.4,3.1,by=0.1)
m2<-m*m

# new variable for tree cover and climate covariate predictions
c<-seq(-4.0,5.0,by=0.1)
c2<-c*c

data.in<-list(spec=nsp,nG=nG,pi=pix,v=v,m=m,m2=m2,c=c,c2=c2,db=db,nsites=nSeg,nVAR=nVAR,y=y,nind=nind,VAR=VAR,dclass=dclass,species=species,
site=site,OBSVAR=OBSVAR,height=height,time=time)
		
inits<-function(){list(N=Nin, mu_a=runif(1,0,1),tau_a=runif(1,0,1),asig=runif(nsp,5,6),mu_s = runif(1, 5, 6),sig_s=runif(1))}
  
params<-c('asig',paste('bsig[', 2:4,']', sep=''),'mu_s', 'sig_s', 'mu_a', 'sig_a', 'mu_b', 'sig_b','Nspec', 'Atot', 'T1', 'T1new', 'Tobs',
	 'Tobsnew', 'Tg1', 'Tg1new', 'r.N', 'alpha', 'beta','Dens','MF.eff','SG.eff','SC.eff','SUN.eff','FFRS.eff','LFP.eff','PAR.eff',
	'Mature.forest','Second.growth','Shade.coffee','Sun.coffee','Fragment','Live.fence','Paramo','Mature','Sec.Growth','Shade.Coff','Sun.Coff',
	'Forest.Frag','Fence','Para','elev.tot','tree.eff','tree.tot','tree.comm','precip.eff','precip.tot','seas.eff','seas.tot',
	'cloud.eff','cloud.tot','Bp.N','Bp.Obs','csig','dsig')

### JAGS model setup 
modelFile='Andean_Community_Model.txt'
parsamples <- run.jags(model=modelFile, monitor=params, data=data.in, inits=inits, n.chains=3, burnin=14000, adapt=1000, thin=5, 
sample=5000, method='rjparallel') 
m1<-summary(parsamples)
m2<-add.summary(parsamples,confidence=c(0.90,0.95)) 




